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ABSTRACT 

We investigate the properties of a closed-form analytic solution recently found by Manko 
et al. (2000b) for the exterior spacetime of rapidly rotating neutron stars. For selected equa- 
tions of state we numerically solve the full Einstein equations to determine the neutron star 
spacetime along constant rest mass sequences. The analytic solution is then matched to the 
numerical solutions by imposing the condition that the quadrupole moment of the numerical 
and analytic spacetimes be the same. For the analytic solution we consider, such a matching 
condition can be satisfied only for very rapidly rotating stars. When solutions to the match- 
ing condition exist, they belong to one of two branches. For one branch the current octupole 
moment of the analytic solution is very close to the current octupole moment of the numeri- 
cal spacetime; the other branch is more similar to the Kerr solution. We present an extensive 
comparison of the radii of innermost stable circular orbits (ISCOs) obtained with a) the ana- 
lytic solution, b) the Kerr metric, c) an analytic series expansion derived by Shibata and Sasaki 
(1998) and d) a highly accurate numerical code. In most cases where a corotating ISCO exists, 
the analytic solution has an accuracy consistently better than the Shibata-Sasaki expansion. 
The numerical code is used for tabulating the mass-quadrupole and current-octupole moments 
for several sequences of constant rest mass. 
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1 INTRODUCTION 

The analytic description of the vacuum spacetime surrounding a 
rapidly rotating neutron star is still an open problem. The analytic 
structure of the spacetime outside a slowly rotating star, and its re- 
lation to the Kerr metric, has been well understood since the sem- 
inal works of Hartle (1968) and Hartle & Thome (1969). On the 
other hand, numerical solutions of the Einstein equations for stars 
rotating up to the mass-shedding limit are now routinely obtained 
with a number of different methods, such as the Komatsu, Eriguchi 
and Hachisu (1989) method (see Stergioulas 2003, for an extensive 
comparison of the different existing numerical methods). These nu- 
merical solutions are indeed useful for modelling astrophysical sys- 
tems, for studying linear perturbations of rapidly rotating relativis- 
tic stars and as initial data for dynamical evolutions of spacetimes 
in numerical relativity (see e.g. Stergioulas & Friedman 1998, Ster- 
gioulas, Kluzniak & Bulik 1999, Stergioulas & Font 2001). 

Despite the availability of numerical solutions, a consistent 
analytic representation of the vacuum metric outside a rapidly ro- 
tating neutron star is desirable for several reasons. In the first place, 
having an analytic form for the metric simplifies the computation of 
the stationary properties of the spacetime. For example, if an accu- 



rate analytic solution were available, geodesies in the neutron star 
exterior could be studied analytically, and one could find closed- 
form expressions for the radii and frequencies of the innermost sta- 
ble circular orbits (ISCOs). In turn, this would simplify the calcu- 
lation of properties of accretion disks, of epicyclic frequencies, of 
accretion luminosities, and so on. 

Furthermore, having an analytic solution could prove use- 
ful to the study of dynamical properties of the spacetime, such 
as gravitational wave emission. One of the unsolved problems in 
gravitational-wave theory is the study of the quasinormal modes 
of rapidly rotating neutron stars. These can be computed either in 
the frequency domain, as an eigenvalue problem, or in the time do- 
main, evolving numerically the (linearized or full) Einstein equa- 
tions and then computing the outgoing radiation. The major tech- 
nical issue in this problem is related to the difficulty of imposing 
outgoing-wave boundary conditions at infinity, since a rapidly ro- 
tating neutron star spacetime is expected to deviate significantly 
from Petrov type D. Having in hand an accurate analytic metric for 
the exterior spacetime one could envisage the possibility of com- 
puting the Weyl scalars in closed form, looking for neutron star 
models which are, in some suitably defined sense, "close to Petrov 
type D" (Baker & Campanelli 2000). If the spacetime is "close 
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enough to type D" one could then apply approximation schemes 
to impose the outgoing-wave boundary conditions. The idea here 
is to improve the presently available methods, which are generally 
based on the use of the Zerilli functions (see e. g. Abrahams et al. 
1992, Allen et al. 1998, Rupright et al. 1998) - i.e., on perturba- 
tions of spherically symmetric vacuum spacetimes. Only recently, 
the Teukolsky formalism for perturbations of Kerr black holes has 
been used for the purpose of wave extraction in the final phase of 
binary black holes mergers (Baker et al. 2002). 

Until the development of a powerful integral equation method, 
devised by Sibgatullin in 1984 (see Sibgatullin 1991 and Manko & 
Sibgatullin 1993 for details), finding analytic solutions to the Ein- 
stein equations for stationary axisymmetric spacetimes was largely 
a matter of guesswork. One typically had to choose some ansatz to 
simplify the mathematical problem of obtaining the solution; then 
one verified a posteriori that the obtained solution had physically 
acceptable properties. In Sibgatullin's method one knows the phys- 
ical characteristics of the solution to be constructed from the very 
beginning, through the choice of the axis expressions of the Ernst 
potentials. 

A complete analytic representation of axisymmetric space- 
times can be obtained in terms of a series expansion whose co- 
efficients are the physical multipole moments (Fodor, Hoenselaers 
& Perjes 1989, Ryan 1995). In principle, this gives an approxima- 
tion to a numerical spacetime that can be made arbitrarily accurate: 
one would need to include a sufficiently large number of multipole 
moments and match them to some given numerical solution. How- 
ever, such a procedure involves a very large number of expansion 
coefficients, which makes it difficult to use for practical purposes. 
Some applications of this idea have already appeared: for example, 
Shibata & Sasaki (1998) derived formulae for the location of the 
ISCO around rapidly rotating neutron stars. 

Quite recently, Manko et al. (2000b) were able to find a new 
asymptotically-flat solution to the Ernst equations for the Einstein- 
Maxwell system. This solution is very interesting because it is 
given in closed form. Furthermore, when two of its parameters 
(i.e., the charge and magnetic moment) are set to zero, the solu- 
tion depends only on three parameters: mass, angular momentum 
and a third parameter b, which is related to the spacetime's physical 
quadrupole moment. With this simplification, the solution reduces 
to a particular three-parameter specialization of the Kinnersley- 
Chitre (1978) solution (a generalization of the Tomimatsu-Sato 
8 = 2 spacetime). Notice however that Kinnersley and Chitre only 
constructed the relevant Ernst potential (they did not provide ex- 
plicit expressions for the corresponding metric functions). Further- 
more, the Kinnersley-Chitre solution is restricted to the subextreme 
case (M 2 > a 2 ). On the other hand, in the solution by Manko et al., 
when electric and magnetic fields are set to zero M and a are al- 
lowed to assume arbitrary real values, because the parameter set in 
their solution is analytically extended. Therefore the Kinnersley- 
Chitre solution is obtained as a particular case of the analytic solu- 
tion in Manko et al. (2000b) when certain restrictions are imposed 
on the parameters of that solution. 

There have been attempts in the literature to fix the free param- 
eters in analytic exterior solutions by matching them to numerical 
solutions. However, different matching conditions were used. For 
example, Sibgatullin & Sunyaev (1998, 2000) fixed the free pa- 
rameters appearing in a different analytic solution using the radii 
of marginally stable circular orbits, or a suitably defined redshift 
parameter at the stellar equator. For their metric, which is different 
from the one we consider here, they found that corrections due to 
the quadrupole moment can accurately reproduce the properties of 



the "exact" exterior spacetime only for several equations of state 
(EOSs), with the exception of EOSs with large phase transitions. 
A simple, closed form expression for the analytic metric used in 
Sibgatullin & Sunyaev (1998, 2000) was given explicitly by Sib- 
gatullin (2002). 

A matching procedure based on the redshift parameter was 
again used by Stute & Camenzind (2002). Our own preference here 
is to avoid matching using local properties and, instead, match the 
solution's mass-quadrupole moment, which is a global property of 
the spacetime. Furthermore, it is well known that deviations from 
the slow-rotation behavior in rapidly rotating stars, due to the stel- 
lar oblateness, are determined mainly by the mass-quadrupole mo- 
ment. The quadrupole moment was also used in matching the ana- 
lytic and numerical solution in Manko et al. (2000a). 

The plan of the paper is as follows. In section 2 we describe 
the procedure to numerically compute the spacetime describing a 
rapidly rotating compact star using the Komatsu-Eriguchi-Hachisu 
(1989) self-consistent field method, as modified by Cook, Shapiro 
and Teukolsky (1994, henceforth CST). In particular, we discuss 
how to implement this method for a numerical evaluation of the 
spacetime's multipole moments. In section 3 we present the an- 
alytic solution recently obtained by Manko et al. (which is only 
valid in the vacuum prevailing outside the rotating neutron star) 
and describe its multipolar structure. In section 4 we describe our 
procedure to match Manko's analytic solution to the numerically 
obtained spacetime, and derive the coordinate transformation re- 
lating the two metrics. Section 5 is devoted to a discussion of the 
tests we used in order to understand "how close" the analytic and 
numerical spacetimes are. As we will discuss in the following, 
there are two possible families of analytic solution for which the 
mass-quadrupole moment of the analytic solution matches to the 
mass-quadrupole moment of the numerical spacetime. The current- 
octupole moment of the first family of solutions is very close to 
the current-octupole moment of the numerical spacetime, while the 
second solution is close to the Kerr spacetime. An examination of 
the metric functions on the equatorial plane and on the rotation axis 
confirms that the first solution is also the one which better approx- 
imates the numerically obtained metric functions. As an indepen- 
dent check, we compute the location of ISCOs in the spacetime sur- 
rounding the rotating star using different approaches. In particular 
we locate ISCOs using the analytic solution, and compare the re- 
sults thus obtained: 1) to the ISCOs found by numerical integration 
of the Einstein equations, and 2) to the analytic formulae for the 
ISCO's obtained by Shibata & Sasaki (1998), truncated at different 
orders of approximation. In most cases where a corotating ISCO 
exists, the analytic solution has an accuracy consistently better than 
the Shibata-Sasaki expansion. Only in some cases the higher-order 
multipoles that are missing in the analytic solution significantly in- 
crease the error in computing the location of the ISCO. Finally, we 
compare our matching procedure to previous work by Manko et al. 
(2000a) and by Stute & Camenzind. The conclusions follow. 



2 NUMERICAL GRAVITATIONAL FIELD OF A 
RAPIDLY ROTATING NEUTRON STAR 

To begin with, in this section we briefly discuss the procedure for 
obtaining highly-accurate numerical solutions for the spacetime of 
rapidly rotating neutron stars and for computing their multipole 
moments. For more details the reader is referred to the review arti- 
cle by Stergioulas (2003). 
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2.1 Numerical determination of the spacetime and 
computation of the multipole moments 

The interior and exterior spacetime of a stationary, axisymmetric 
star is described by a metric in the following form: 

2 - -e 2v dt 2 +fiV 2 Vsin 2 9(d<j)-C0A) 2 + e 2a (dr 2 + r 2 dQ 2 ), 



ds 



where v, B, a and CO are four metric functions to be determined by 
solving four field equations. In the numerical method of Komatsu 
et al. (1989, henceforth KEH) one defines two auxiliary functions 
p, y through the relations v = (y + p)/2 and B = e 1 . Then, three 
out of the four field equations are written in the following integral 
forms 

PM=-L e~ V2 

n=0 

U~dr'j\jr a fl{r,r')P 2n {J)S- p {riJ)} 

Plnin), (2) 



For a configuration that is stationary, axisymmetric, sym- 
metric with respect to reflections in the equatorial plane and 
asymptotically flat, the spacetime can be characterized by two 
sets of scalar multipole moments: the even- valued mass mo- 
ments (Mo, M 2 , M4...) and the odd- valued current moments 
(1) (Si , 53 , 55 . . .)• Ryan (1997) presented a method for extracting the 
multipole moments from the asymptotic form of the metric func- 
tions. The lowest-order appearance of each moment in terms of a 
power series in 1 jr is determined by the expansions 
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and Sp, Sy and S a are lengthy source terms, whose expressions can 
be found in KEH. In the equations above, /j = cos 6, while P n (ju) 
denotes the Legendre polynomials and P™(/j) the associated Leg- 
endre functions. The metric function a is determined by an ordinary 
differential equation. 

We compute numerical equilibrium models using the code by 
Stergioulas and Friedman (1995) (see Nozawa et al. 1998 and Ster- 
gioulas 2003 for extensive accuracy tests). The numerical code uses 
the CST formulation, in which the KEH equations are written in 
terms of a compactified coordinate s defined through the relation 



s 



(5) 



where r e is the (coordinate) radius of the stellar equator. This al- 
lows the computation of the whole exterior spacetime out to infin- 
ity, which is important in detailed comparisons of the numerical 
metric to the analytic metric. 
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By comparison with equations (2) and (4) one finds that 

M 2 „ = \ j~dr' jT 1 d t lr' 2n+2 P 2n ( f l)S- p {r',iJ), (8) 



(9) 



and 

Sm-x = ±£d, / £d f // 2 * 2 smtfPl,_ l (jl)S a (T', f /). 

Thus, in general, any model of a rapidly rotating neutron star has an 
infinite number of mass- and current-multipole moments. In order 
to match an analytic exterior metric to a numerically-computed in- 
terior metric and to check the accuracy of the matching procedure, 
we computed the mass-quadrupole moment M 2 and the current- 
octupole moment S3 . 

An alternative, asymptotic method for evaluating the multi- 
pole moments was introduced by Laarakkers & Poisson (1997). We 
also used their method in order to cross-check the results obtained 
from the integral relations (8) and (9). The idea, in this case, is to 
evaluate numerically the coefficient of P 2n (M) in the general expres- 
sion for (2) - or analogously, the coefficient of P\ n _ l in (4) - at the 
outermost grid points (i.e., as °o), and multiply the result by 
the appropriate factor (containing powers of r) that can be obtained 
from equations (6) and (7). We have checked that the two methods 
typically agree to better than one part in 10 3 . 



2.2 Equilibrium sequences 

The equilibrium solutions for a given EOS form a two-parameter 
family. In particular, stable equilibrium solutions are bounded by 
four limit sequences. These limits are shown in Figure 1, which 
displays the gravitational mass M vs. the central energy density e c 
for one the EOSs derived by Akmal, Pandharipande and Ravenhall 
(1998, henceforth APR). As an illustrative example we consider 
the APR EOS which does not include boost interactions (we refer 
to the original paper for details). The qualitative picture does not 
change when we consider other EOSs (see eg. CST, where plots 
are presented for a representative sample of EOSs). The solid line 
is the static limit - that is, the sequence of nonrotating solutions 
to the standard Tolman-Oppenheimer-Volkoff equations. The long- 
dashed line is the mass-shedding (Kepler) limit, which is deter- 
mined by the condition that the centrifugal force exactly balances 
the gravitational attraction at the stellar equator, in which case a 
fluid element on the equator has the same angular velocity as a free 
particle in a Keplerian orbit at the same location. Both sequences 
terminate at high central density at the stability limit, where equi- 
librium solutions are marginally stable to axisymmetric perturba- 
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Figure 1. Limit sequences for EOS APR: the solid line corresponds to the 
nonrotating limit; the long-dashed line corresponds to the mass-shedding 
(Kepler) limit; the dotted line is the axisymmetric instability limit. The 
nearly horizontal lines are sequences of constant rest mass. From bottom to 
top: the dashed line corresponds to a star of gravitational mass M = 1 AM e 
in the non rotating limit; the dash-dot line is a maximum-mass normal se- 
quence; and finally, the dash-dot-dot line is a selected supramassive se- 
quence. 



tions; and they terminate at low central densities at the low-mass 
limit, below which a neutron star cannot form (not shown in Figure 
1). 

Within the class of stable equilibrium solutions, CST pointed 
out the significance of constant rest-mass sequences, called evo- 
lutionary sequences, since an isolated neutron star, slowly losing 
energy and angular momentum via some dissipative process (e.g. 
electromagnetic emission or gravitational-wave radiation), must 
evolve conserving the total baryon number, and hence its rest mass 
Mb- An accreting neutron star in a binary system will not evolve 
along a constant rest-mass sequence: the actual sequence depends 
on several parameters, such as the magnetic field, accretion rate etc. 
Nevertheless, the constant rest-mass sequences in CST have been 
used in the past in evaluating the accuracy of analytic exterior so- 
lutions and we will also use them here solely for the same reason. 
We compute three constant rest-mass sequences for each EOS: 

• the sequence corresponding to a canonical neutron star having 
gravitational mass M = 1 AM® in the non rotating limit, 

• the sequence terminating at the maximum-mass model in the 
non rotating limit (maximum-mass normal sequence). 

• a supramassive sequence, i.e., a sequence which does not ter- 
minate at a nonrotating model. 

We include the following set of EOSs. For comparison with 
CST we include EOSs A, AU, FPS and L. We refer to their pa- 
per for an extensive discussion of each EOS. We supplement the 
set of EOSs considered by CST with a relatively new model: the 
model derived by Akmal, Pandharipande and Ravenhall (1998) 
from Hamiltonian many-body theories of nuclear matter, includ- 
ing boost corrections in the Hamiltonian (henceforth, we will refer 
to this model as APR-b, where "b" stands for "boosted"). 

In Tables (1-5) we give numerical results for the structure 
properties of the models we have computed. All models have been 
computed using a resolution of (301 angular points) x (601 ra- 
dial points), corresponding to a typical accuracy of at least one part 
in 10 3 in all quantities. Each Table corresponds to a constant rest 
mass sequence, and lists: the total central energy density e e in units 
of 10 15 g citT 3 ; the angular velocity Q. in units of 10 3 s _1 ; the mo- 



ment of inertia / in units of 10 45 g cm 2 (for rotating models only); 
the gravitational mass M in solar masses; the ratio of rotational 
kinetic energy to gravitational binding energy T/W; the equatorial 
circumferential radius of the star R e in km and the height (in km) of 
corotating (h+) and counterrotating (h-) ISCOs from the surface of 
the star (if an ISCO does not exist, the corresponding entry is omit- 
ted). The height of an ISCO is defined as the difference between the 
circumferential radius at the ISCO and the circumferential equato- 
rial radius of the star. The next three columns give the first few 
physical multipoles in geometrized units of c = G = 1 : namely, we 
list the mass quadrupole moment M2 = Q in km 3 , the angular mo- 
mentum S\ = J in km 2 , and the current octupole moment 53 in 
km 4 . We have checked our code by reproducing the quadrupole 
moments computed by Laarakkers & Poisson and found excellent 
agreement. The accuracy in computing 53 was checked by compar- 
ing the integral form to the asymptotic form mentioned in Section 
2.1, finding good agreement. This shows that using the compacti- 
fied coordinate introduced in CST allows a very accurate numerical 
determination of relatively high-order multipoles. 

The last column gives the value of the "quadrupole" parameter 
b (in km) for which the analytic solution provides a good approxi- 
mation of the numerical spacetime. When matching the quadrupole 
moment of the numerical and analytic spacetimes is not possible, 
the corresponding entry is omitted. More details on the procedure 
we followed to obtain the values listed in this column will be given 
in section 4. 



3 ANALYTIC GRAVITATIONAL FIELD OF A RAPIDLY 
ROTATING NEUTRON STAR 

In this section we summarize the properties of the vacuum analytic 
solution obtained by Manko et al. (2000b). We will concentrate 
in particular on the multipolar structure of the solution, since our 
ultimate purpose will be to reproduce accurately the first few mul- 
tipoles of the numerical spacetimes we discussed in the previous 
section. 



3.1 The solution by Manko et al. 

In the vacuum region surrounding a stationary and axisymmetric 
star, the spacetime only depends on three metric functions (while 
four metric functions are needed for the interior). The most general 
form of the metric (Papapetrou 1953) is 

ds 2 = -f(dt-wd§) 2 +f- x {e 2 ~<(dp 2 +dz 2 ) + p 2 d<$> 2 } . (10) 

Here /, w and J are functions of the quasi-cylindrical Weyl-Lewis- 
Papapetrou coordinates (p, z). Starting from this metric one can 
write down the vacuum Einstein-Maxwell equations as two equa- 
tions for two complex potentials •£ and <I>, following a procedure 
due to Ernst (1968). The equations are: 

(tfe{£} + |*| 2 )V 2 £ = (VE +2**V*)-VE (11) 
(^e{E} + |*| 2 )V 2 * = (VE + 24>*V<t>)V<I> (12) 

Once the potentials are known, the metric can be reconstructed. 
Sibgatullin (1991) devised a powerful procedure for reducing the 
solution of the Ernst equations to simple integral equations. Basi- 
cally, one starts with a choice for the values of the Ernst potentials 
on the symmetry axis 

e(z) = £(j5 = 0,z), /(2) = *(p = 0,f), (13) 
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solves two complex-valued integral equations, and checks that the 
obtained solution satisfies the expression of the Ernst potentials in 
terms of physical multipoles: 



1-5 
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The real parts of m n are the mass multipoles, the imaginary parts 
of m n are the current multipoles, the real parts of q„ are the electric 
multipoles and the imaginary parts of q„ are the magnetic multi- 
poles. 

After more than ten years of work in the field, Manko et 
al. (2000b) were finally able to find a vacuum solution involving 
five parameters (mass, angular momentum, charge, magnetic dipole 
moment and mass quadrupole moment) which can be expressed in 
terms of relatively simple rational functions. We are particularly in- 
terested in solutions having no charge or magnetic dipole moment. 
If we denote by M the gravitational mass of the star, by a the spe- 
cific angular momentum {a = J/M), and introduce a parameter b 
which can be related to the mass quadrupole moment, their choice 
for the axis values of the Ernst potentials is: 
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(19) 



To be able to write the metric in rational form, one must introduce 
generalized spheroidal coordinates 

r+ + r_ 



r+ — r- 

A= 2k ' y= 2k ' 
where r± = \/p 2 + {z±k) 2 and 

k= y/d + 8. 



(20) 



(21) 



The inverse transformation between the two sets of coordinates is 
p = fc(l-y 2 ) 1 / 2 (x 2 -l) 1 / 2 , z = kxy. (22) 
The metric is then written as 



ds z 



with 
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(23) 
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and 



D = {A(k 2 x 2 - Sy 2 ) 2 + 2kmx[2k 2 (x 2 - 1) 

+ (28 + ab-b 2 )(l-y 2 )} 

+ (a-b)[(a-b)(d-8)- m 2 b] (y 4 - 1) - Ad 2 } 2 

+ Ay 2 {2k 2 (x 2 - 1 ) [kx(a - b) - mb] - 2mb8( \-y 2 ) 



+ [(a-b)(d-8)-m 2 b](2kx + m)(l-y 2 )} 2 

= {A[k 2 (x 2 -l)+8(l-y 2 )} 2 
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- I6k 2 (x 2 -i)(l-y 2 ){(a-b){k 2 {x 2 -y 2 )+ 28y 2 ] 



(25) 
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F = Sk 2 {x 2 -l){(a-b)[k 2 {x 2 -y 2 )+28y 2 ]+y 2 m 2 b} 

x {kmx[(2kx + m) 2 - 2y 2 (28 + ab - b 2 ) - a 2 + b 2 ] 

- 2y 2 (A8d-m 2 b 2 )} + {A[k 2 {x 2 -l)+8(l-y 2 )] 2 
+ (a-b)[(a-b){d-8)-m 2 b](l-y 2 ) 2 } 

x (A(2kmbx + 2m 2 b)[k 2 (x 2 -l) + 8(l-y 2 )} 

+ (\-y 2 ){(a-b)(m 2 b 2 -AM) 

- (Akmx + 2m 2 )[(a-b)(d-8)-m 2 b]}). 



(27) 



In order for the solution to satisfy the requirements of axisymme- 
try, stationarity and reflection-symmetry in the equatorial plane, all 
three parameters M, a and b must be real. 

3.2 Multipolar structure of the analytic solution 

Here we examine the multipolar structure of the analytic solution 
by Manko et al. for rotating and nonrotating solutions. The only 
nonvanishing multipole moments of the solution are the gravita- 
tional mass Re{mo} = M, the quadrupole moment Re{ni2} = Q, 
the angular momentum Im{m\ }=J = aM and the current octupole 
lm{m-i\ = 53. The quadrupole moment and the current octupole 
moment are given in terms of the three parameters M, a and b as 

Q = -M(d-8-ab + a 2 ), (28) 

and 

S 3 = -M<[a 3 -2a 2 fe + a[ft 2 + 2(d-8)] -ft(d-S)}. (29) 

However, since a and b are independent parameters, setting a equal 
to zero does not automatically imply a vanishing Q and 53 , as would 
be the case for a realistic solution of a nonrotating perfect fluid star. 
Instead, the nonrotating solution (a = 0) has a quadrupole moment 
equal to 

u2\ 2 



Q(a = 0) 



M (M 



A (M 2 -b 2 ) ' 
and a current octupole moment equal to 

S 3 (a = 0) = -Z?2(a = 0). 



(30) 



(3D 



It is obvious that there exists no real value of the parameter b for 
which the quadrupole moment vanishes for a nonrotating star. For 
\b\ <M, the solution is oblate (Q < 0) with a minimum quadrupole 
deformation obtained for b = 



lelmin(« = 0)=M 3 /4. 



(32) 



At b = ±M, the nonrotating multipole moments Q and 53 diverge, 
while for \b\ > M, the nonrotating solution is prolate (Q > 0) with 
a minimum quadrupole deformation of 

emin(« = 0)=2M 3 , (33) 
at b = ±V3M. 

Obviously, the analytic solution by Manko et al. does not re- 
duce continuously to the Schwarzschild solution as the rotation 
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Figure 2. The matching condition Q - Qn = as a function of b for the 
fastest rotating FPS model in the maximum-mass sequence (last row of 
M B = 2.1O5M.0 sequence in Table 3). There are two possible solutions: a 
solution for which (typically, as in the case shown) b = b- < 0, and a pos- 
itive solution for which b = b + > 0. The "negative" solution b = b- is the 
one which is relevant for rapidly rotating neutron stars, as shown in section 
5. 

vanishes. It can only reduce to other forms of nonrotating vac- 
uum solutions (e.g. the well-known Weyl solutions) that could be 
matched to other interior solutions, such as nonrotating stars with 
non-isotropic stresses, inducing nonvanishing quadrupole deforma- 
tions. Nevertheless, as we will show next, the analytic solution can 
approximately describe a rapidly rotating fluid star, when the ro- 
tation rate is large enough, so that the quadrupole deformation in- 
duced by the rotation roughly exceeds the minimum nonvanishing 
oblate quadrupole deformation of the solution in the absence of ro- 
tation, i.e. roughly when 

\Q\>M 3 /4. (34) 

Since the quadrupole moment is roughly proportional to a 2 M, one 
expects that the analytic solution could be relevant for rotation rates 
of roughly j > 0.5, where j = J/M 2 is a dimensionless measure of 
the angular momentum of the star. We will confirm this expectation 
by direct comparisons with numerical solutions in the next section. 

It is interesting that the Kerr solution can still be obtained from 
the analytic solution, if one accepts the following imaginary form 
for the parameter b 

b = i\/M 2 -a 2 , (35) 

with a < M. In this case, one recovers the correct expressions Q = 
-a 2 M and S 3 = -c?M. 



4 MATCHING THE INTERIOR AND EXTERIOR 
SOLUTIONS 

The three parameters of the analytic solution (M,a and b) can be 
set at will. However, only certain combinations of values can corre- 
spond to specific models of rapidly rotating neutron stars. Since the 
solution has four nonvanishing multipole moments, but only three 
free parameters, one can at most match three multipole moments 
of any given numerical solution. The fourth multipole moment will 
then be determined by the analytic solution, and its relative differ- 
ence with the (known) numerical value will be a measure of the 
accuracy of the analytic solution. The four multipole moments are 
not equally important for specifying a solution, 53 being the least 




-6.0 -4.0 -2.0 0.0 2.0 4.0 6.0 
b 



Figure 3. The matching condition Q - Qn = as a function of b for a slowly 
rotating model in the maximum-mass sequence for EOS FPS (fourth row of 
M B = 2.105M Q sequence in Table 3). No real solution exists. 

important, even for the most rapidly rotating models. Therefore we 
choose to match the analytic exterior solution to known numerical 
solutions by matching the gravitational mass M, the specific an- 
gular momentum a and the mass-quadrupole moment Q. One then 
hopes that the resulting analytic solution will yield a value for the 
current-octupole moment 53 that is close to the corresponding value 
in the numerical model. As we will show, there exists a branch of 
solutions for which this is indeed the case. Manko et al. (2000a) 
also used the quadrupole moment to match numerical and analytic 
solutions, but their examples correspond to the other branch of so- 
lutions, for which the analytic value of S3 does not agree well with 
the numerical value. 

For a given model of a rapidly rotating neutron star, we first 
construct a highly accurate numerical solution, as described in sec- 
tion 2. In the analytic solution (23), we set M and a to be equal to 
the obtained numerical values. The remaining parameter b is then 
determined by solving the equation 

Q-Qn=0, (36) 

where <2n is the value of the quadrupole moment obtained by the 
numerical code. A plot of Q — <2n as a function of the parameter b, 
for the most rapidly rotating model of the maximum-mass sequence 
for EOS FPS, is shown in Fig. 2. Two possible real solutions for b 
exist: a solution that is usually negative, b-, and a solution that 
is always positive, b + . Thus, for each set of physical parameters 
M,a and Q, there exists two different branches of solutions, with 
parameters (M, a, b- ) and (M, a, b+), respectively. In the remainder 
of the paper, we will refer to these two different branches as the 
negative solution (-) and the positive solution (+). As we will show 
next, these two branches correspond to very different spacetimes. 

In the previous section, we estimated that the analytic solution 
should be relevant for rapidly rotating neutron stars only for values 
of j roughly larger than 0.5. Fig. 3 shows a more slowly rotating 
model than the model shown in Fig. 2, along the same evolutionary 
sequence. It is obvious that no solution to equation (36) exists for 
any real value of the parameter b. Along each sequence there is 
a critical rotation rate above which one can match the numerical 
interior solution to the analytic exterior solution. In Tables (1-5) 
we list all computed physical properties for the selected sequences. 
The last column lists the parameter b = b- of the negative branch 
of the analytic solution (when it exists). This is the relevant branch 
for rapidly rotating neutron stars, as we will show in section 5. In 



Approximate Matching of Analytic and Numerical Solutions for Rapidly Rotating Neutron Stars 7 



Tables (1-5) some models appear having b = b- > 0. These models 
do not belong to the positive branch b + . They are instead models 
which are very close to the critical value of the rotation parameter, 
j = j cr i t : in these particular cases, both solutions to equation (36) 
can happen to be positive. However, in general (as long as a model 
is rotating somewhat above the critical rate) the negative branch has 
< 0. 

Typical values of j cra above which the analytic solution ex- 
ists are listed in Table 6 for a subset of the considered EOSs. For 
smaller masses, j cr i t is usually smaller: therefore, for "canonical" 
neutron stars (having mass M ~ 1 AMq in the non rotating limit) 
the analytic solution is valid over a wider range of j. In terms of 
the angular velocity at the mass-shedding limit for uniformly rotat- 
ing stars, the critical rotation rates are given in the right column of 
Table 6. The critical rotation rate £l C rit /^-Kepler ranges from ~ 0.4 
to ~ 0.7 for the M = 1 .4M Q sequence, with the lower ratio corre- 
sponding to the stiff est EOS. For the maximum-mass sequence the 
ratio is ~ 0.9, nearly independent of the EOS. In conclusion, the an- 
alytic exterior solution can be useful for studying rapidly rotating 
neutron stars. The exterior gravitational field of massive neutron 
stars created in binary neutron star mergers, supported temporarily 
by differential rotation against collapse, could also be described, to 
some accuracy, by the analytic solution (the accuracy will depend 
on how significant the higher multipole moments are in the case of 
strong differential rotation). If the EOS is very stiff, such as EOS 
L, then the analytic solution is also valid for for description of ac- 
creting neutron stars in Low-Mass-X-Ray binaries (LMXB), with 
rotational periods of a few milliseconds. 

4.1 Coordinate transformations between vacuum and 
non-vacuum metrics 

Before presenting specific tests of the accuracy of the analytic solu- 
tion, we need to describe the coordinate transformation that relates 
the interior metric (1) to the exterior vacuum metric (10) (see Is- 
lam 1985) . For the interior metric (1), we define the cylindrical 
coordinates 
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Figure 4. The current octupole moment 53 as a function of j for the numer- 
ical solution, for the Kerr metric, and for the negative and positive branches 
of the analytic solution. For illustrative purposes we have chosen the FPS 
EOS, and fixed our attention on the sequence having maximum mass in the 
non rotating limit. The negative branch reproduces with excellent accuracy 
the numerical behavior of the current octupole, so it is the branch appropri- 
ate for describing the exterior spacetime of rapidly rotating neutron stars. 
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the metric in the exterior takes the desired form (10). Since the 
transformation (41,42) for the coordinate z cannot, in general, be 
solved analytically, one can relate a solution for the interior metric 
(1) to the exterior metric (10) only through numerical integration. 



O5 = rsin0, z=rcos9. 

In vacuum, Einstein's field equations imply that 

0. 



d 2 (m) d 2 (m) 



dm 2 dz 2 
One can therefore define a new coordinate 

p = 03S, 



(37) 



(38) 



(39) 



satisfying the two-dimensional Laplace equation (38), and a second 
coordinate 



z = z(0J,z), 

satisfying the Cauchy-Riemann conditions 



3z 
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dz 
dz 



dz 
3p = 
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The coordinate z is obtained by integration of the above Cauchy- 
Riemann conditions, requiring that z = in the equatorial plane (at 
z = 0). It is easy to show that 

-1 



d® 2 +dz 2 ■ 



dp 1 



-dz 2 



(43) 



5 TESTS OF THE ACCURACY OF THE ANALYTIC 
SOLUTION 

5.1 The current-octupole moment 

The current-octupole moment 53 , like the quadrupole moment Q = 
Mi, is a function of a, b and M. Once we have fixed b by matching 
the quadrupole moment to the numerical spacetime through equa- 
tion (36), there are no more free parameters to be specified; the 
current-octupole 53 can be computed using equation (29), and then 
compared to the value of 53 computed for the numerical metric. 
Therefore, 53 serves as a good error indicator for the accuracy of 
the solution. In fact, the value of 53 obtained analytically for the 
two branches of solutions, b + and b- , can be used to distinguish 
which of the solutions is more relevant for rapidly rotating neutron 
stars. Fig. 4 displays 53, for the two branches of the analytic so- 
lution, along with the value of 53 for the numerical solution and 
for the Kerr solution, for the evolutionary sequence corresponding 
to EOS FPS and having maximum mass in the nonrotating limit 
(see Table 3). The error for the (-) solution is very small, at most 
of the order of 3 %. On the other hand the error for the (+) solution 
is quite large (up to 56 % for the fastest rotating model). In this 
case the solution is closer to the Kerr value than to the value cor- 
responding to numerical models of rapidly rotating neutron stars. 
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Figure 5. Relative error in the g„ -component of the analytic metric and of 
the Kerr metric in the equatorial plane, when compared to the numerical 
solution for a rapidly rotating star. 

Therefore, in the remainder of this paper we will only use the (-) 
branch of solutions to the matching condition (36). 

In Table 7 we display the relative error AS in the analytic 
value of 53 (when compared to the numerical solution) for the crit- 
ically rotating and maximally rotating models of all evolutionary 
sequences considered in this paper. For most sequences, the relative 
error in 53 can be as large as 12% for the critically rotating models, 
reducing to a few percent only for the models at the mass-shedding 
limit. Typically, the largest errors appear for the 1 AMq sequences, 
while the sequences that terminate at the maximum mass static 
model have the smallest errors. In most cases the error is larger 
for slower rotating models. This shows that for those models 53 is 
still influenced by its nonzero value in the non rotating case (for 
the analytic solution, the octupole moment 53 , like the quadrupole 
moment Q, does not vanish for a = and b 7^ 0). For more rapidly 
rotating models this influence diminishes and 53 becomes almost 
entirely of rotational origin, agreeing better with the numerical so- 
lution. Comparing the various EOSs, one sees that the error in 53 
for soft EOSs, such as EOS A, is smaller than the corresponding 
error for very stiff EOSs, such as EOS L. The critical 1.4M S model 
for EOS L shows an unusually large relative error of 45% in 53. 
This, again, is related to the compactness of the various models and 
to the value of the multipole moments for a = 0. 

5.2 Direct comparison of metric components 

As a second test of the accuracy of the analytic solution for rapidly 
rotating neutron stars, we performed a direct comparison of spe- 
cific metric components for several representative models, using 
all EOSs in our sample. Here, we focus on the most rapidly rotat- 
ing model of the maximum mass sequence with EOS FPS, since 
the other cases we examined showed similar behavior. 

For this model we computed the metric components g tt , g t § 
and gAA on the equatorial plane and along the symmetry axis using 
the analytic metric and the Kerr metric. Then we compared the rel- 
ative error of both metrics with respect to the corresponding com- 
ponents of the numerical metric. 

Fig. 5 shows the relative error of the g„-component of the an- 
alytic metric and of the Kerr metric in the equatorial plane. For 
the analytic metric the error is only 0.3% at the surface of the star 
(located at p = 10.6), and decreases monotonically with increas- 
ing distance, becoming of order 10~ 6 near infinity. In comparison, 



Figure 6. Same as Fig. 5, but along the axis of symmetry. 

the relative difference between the Kerr metric and the numerical 
metric is 1.3% at the equator (i.e., four times larger than the er- 
ror in the analytic metric). The relative difference between the Kerr 
metric and the numerical metric also decreases with increasing dis- 
tance, as expected, and for distances larger than about 200 equato- 
rial radii, the difference between the analytic solution and the Kerr 
solution is negligible. In other words, at such a distance the error in 
the analytic solution is dominated by the Kerr contribution at first 
order in the rotational parameter, while the effects of the higher- 
order multipole moments Q and 53 have become unimportant. The 
corresponding figure for the relative error in on the equatorial 
plane is nearly identical to Fig. 5 for g tt . When we consider g,^ in 
the equatorial plane, the relative error at the surface is 1.3% for the 
analytic metric and 5.3% for the Kerr metric. This larger error for 
g,§ should be expected: this metric component vanishes in the non- 
rotating limit, so it is more sensitive to contributions by the higher- 
order multipole moments Q and 53 than the metric components g t , 
andg^. 

In order to compare the metric components on the symmetry 
axis, we first need to integrate the Cauchy-Riemann conditions (41) 
and obtain the coordinate z in terms of the coordinate z. This can 
be done easily, once the numerical solution for the metric function 
B is obtained. Fig. 6 shows the relative error in g tt for the analytic 
solution and the Kerr solution along the symmetry axis. The loca- 
tion of the surface (as determined from the numerical solution) is at 
z = 6.05. At the surface, the relative error for the analytic solution 
is 7%, while it is 15% for the Kerr metric. Thus, we see that the 
effect of a large quadrupole moment Q shows up predominantly 
in the metric components along the symmetry axis, while in the 
equatorial plane this effect is very small. The reason for this differ- 
ence is that a rapidly rotating star becomes very oblate, so that the 
stellar surface on the symmetry axis is located deeper in the gravi- 
tational potential well than the surface in the equatorial plane. The 
specific example shown in the above figures has polar to equato- 
rial axes ratio of 0.6, thus the equatorial radius is roughly twice as 
large as the polar radius. The analytic value of g n on the surface in 
the equatorial plane is —0.59, while it is —0.44 on the surface on 
the symmetry axis (the asymptotic value at large distances is — 1). 
Gravity is stronger on the polar surface, and this justifies a larger 
relative error in g tt there. At about 3 polar radii, the relative error in 
gu along the symmetry axis decreases to the 1% level for both the 
analytic and Kerr solutions. 

The above direct comparison of metric components shows that 
the analytic metric is a good approximation to the numerical one 
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(or, at least, a much better approximation than the Kerr metric) in 
the equatorial plane, where one expects particle orbits to be as- 
trophysically more relevant. For gravitational-wave extraction in 
numerical relativity, the larger inaccuracies near the polar surface 
could influence the waveforms. In order to minimize this effect, the 
extraction should be done as far as possible from the surface of the 
star. In any case, the analytic metric is everywhere more accurate 
than the Kerr metric. Thus a perturbative wave-extraction scheme, 
built with the analytic metric as a background, should yield more 
accurate waveforms than those obtained with techniques available 
at present (which are based on a perturbative extraction of wave- 
forms around a Schwarzschild or Kerr background). 

5.3 Innermost stable circular orbits 

It is well known that not all orbits around relativistic stars are sta- 
ble. For nonrotating stars, the ISCO is located at a circumferential 
radius of Rjsco = 6M . Depending on the EOS and the mass of 
the star, the ISCO can be located outside the stellar surface. Rota- 
tion introduces a preferred direction in the (]) coordinate, so ISCOs 
around a rotating star belong to two distinct families: a corotating 
and a counterrotating one. For moderate rotation rates, the effect 
of rotation is to shorten the distance between the surface and the 
corotating ISCO. For rapid rotation, the large quadrupole moment 
of the star reverses this trend (notice that even in some Newtonian 
stellar models, large higher-order multipole moments can introduce 
an ISCO (see, Zdunik & Gourgoulhon, 2001; Amsterdamski et al. 
2002). The counterrotating ISCO radius normally increases with 
rotation. Detailed computations of ISCOs for a large number of 
models and EOSs are presented in CST; we also refer the reader to 
that paper for the equations defining the ISCOs that were used in 
our numerical computations. 

Testing the accuracy of the analytic solution in computing the 
properties of ISCOs is important, as ISCOs are related to several as- 
trophysical properties of rapidly rotating neutron stars in LMXBs. 
An accretion disk cannot extend to radii located within the ISCO, 
and this sets an upper limit to the Keplerian frequency of parti- 
cles orbiting a star. This idea could be used, e.g., in determining 
whether compact stars in LMXBs are composed of strange mat- 
ter (Stergioulas, Kluzniak & Bulik 1999, Gondek et al. 2001). In 
addition, the location of the ISCO could play a role in the mech- 
anism producing the kHz quasi-periodic oscillations observed in 
many LMXBs (van der Klis 2000): see, e.g., Kluzniak et al. (2003). 

A circular orbit in the equatorial plane is one for which 05 = 
const., and hence p = const. The equation for geodesic motion 
along the radial coordinate p reads 

where E and L are the conserved energy and angular momen- 
tum per unit mass, determined by the conditions V = dV /dp = 0. 
Geodesies become unstable when d 2 V/dp 2 = 0, or 

(wV'/ 5 p(2/-/'p) + w' 2 / 4 [2/ 2 + (-/' 2 +f"f)p 2 ] 

W/ 2 ^w' 2 / 4 + /'p(2/-/'p) 

[2/ 2 + 2/' 2 p 2 - /p(4/' + f"p)] + p(2/ - f'p) 
{3/7 2 -4/' 2 /p + /' 3 p 2 +/ 2 [/"p 

-wVvW 4 +/'P(2/-/'P)]})/ 



(/ 2 p 2 {w' 2 / 4 +3/7p-/' 2 p 2 

-/ 2 [2 + w' v W 4 +/'p(2/-/'p)]}) =0 (48) 

for corotating orbits (cf. Stute & Camenzind, 2002), where ' in- 
dicates a partial derivative with respect to p. For counter-rotating 
orbits, one can simply use the above equation and change the sign 
of the star's angular momentum. 

Shibata & Sasaki have used a more general representation of 
axisymmetric vacuum solutions (in the form of a series expansion 
that is completely determined by the physical multipole moments 
of the spacetime: see Fodor, Hoenselaers and Perjes 1989 and Ryan 
1995) and derived an approximate analytic formula for the location 
of the ISCO. Their formula depends on the stellar mass, angular 
momentum, mass quadrupole, current octupole and mass 2 4 -pole 
moments. Including all terms up to order 0(4) in the rotation pa- 
rameter, they find the following equation for the circumferential 
radius of the corotating ISCO: 

Rjsco = 6M( 1- 0.54433 j- 0.22619/ +0. 17989£ 2 

- 0.23002j 3 +0.26296j2 2 -0.05317 g3 

- 0.29693 / + 0.44546/ Q 2 - 0.0624902 

+ 0. 01544^4-0. lBlO;"^). (49) 

In the previous expression we have introduced dimensionless pa- 
rameters Q 2 = -Q/M 3 , q 3 = -S3/M 4 and Q A = M4/M 5 . In the 
case of the Kerr metric, the approximate expression for the loca- 
tion of the corotating ISCO up to order O(f) is (see, e.g., Shapiro 
& Teukolsky 1983) 

R fsco = 6M( 1-0.54433; -0.04630/ (50) 

- 0.02016/ -O.OIIIO/). 

The location of the counter-rotating ISCO is obtained from the 
above formulae by reversing the sign of j and 53 . 

For illustrative purposes we again focus on the three sequences 
of EOS FPS (Table 3). We have carried out the calculation for other 
EOSs as well. Although there are quantitative differences between 
the various models, the qualitative behavior and the relative accu- 
racy between the numerical and analytic solutions remain similar 
to those shown here. 

For each sequence we find the relative error in computing 
corotating and counter-rotating ISCO radii with respect to the nu- 
merical solution. We perform this comparison for the analytic solu- 
tion obtained through our matching procedure (when a solution to 
the matching condition exists), for the Shibata-Sasaki formula (49) 
and for the Kerr formula (50). In the case of the Shibata-Sasaki 
formula we do not explicitly compute the moment M4 from the nu- 
merical solution, but we follow the same approximation adopted by 
Shibata & Sasaki. Namely, we set Q4 = OI4Q 2 , where 0C4 is expected 
to take values ranging between and 2. Again, following Shibata 
& Sasaki, we normally set 0C4 = 1 (unless otherwise noted). 

Fig. 7 shows the relative errors in computing the ISCOs for the 
sequence having M = 1 .4M Q in the non rotating limit. In this and in 
the following Figures, negative values of j correspond to counter- 
rotating orbits, while positive values of j correspond to corotating 
orbits. In the corotating case the ISCO disappears at slow rotation 
rates, even before the analytic solution becomes valid; therefore, 
for this sequence, we can only compare the accuracy in finding 
counter-rotating ISCOs. For the Kerr solution the error increases 
monotonically with | j\ , reaching 1 1 % for the fastest rotating model. 
On the other hand, the error for the Shibata-Sasaki formula with 
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0C4 = 1 is only 2% for the fastest rotating model. For this sequence 
the analytic solution is initially close to the Shibata-Sasaki formula, 
but then shows a rather large error, that becomes 10% for the fastest 
rotating model. The explanation for this behavior is that at very 
large rotation rates the inclusion of the multipole moment M4 is 
important, but this multipole moment is absent in the analytic solu- 
tion. When we omit M4 in the Shibata-Sasaki formula (while still 
keeping all other mixed terms up to order 0(4)) by setting (X4 = 0, 
we obtain an error which is much closer to the error made using 
the analytic solution. On the other hand, including only terms up 
to 0(3) in the Shibata-Sasaki formula gives a much smaller error, 
comparable to (or better than) the error of the formula when all or- 
ders up to 0(4) with 0C4 = 1 are included. What this comparison 
underlines is the importance of being consistent up to a certain or- 
der in the rotation parameter. The Shibata-Sasaki formula has small 
error when used consistently up to 0(3) or up to 0(4), but a large 
error when only a few mixed terms up to 0(4) are included. The 
error of the Shibata-Sasaki formula to order 0(4) should improve, 
if one would include the precise values for M4, instead of the crude 
estimate of (X4 = 1 . The analytic solution suffers from the inconsis- 
tency that while the M4 moment vanishes, it is still an exact ana- 
lytic solution. This means that mixed terms containing j, Q and 53 
up to order 0(4) are present. It follows that, for the counterrotating 
ISCOs in Fig. 7, the analytic solution is not as accurate as a con- 
sistent application of the Shibata-Sasaki formula. Notice that the 
non-monotonic increase in the error for the counter-rotating ISCO 
with the Shibata-Sasaki formula at large rotation rates is a con- 
sequence of the moment M4 becoming important near the mass- 
shedding limit: for sequences of larger mass, as the ones we exam- 
ine next, M4 appears to be much less important. 

The comparison of the error in computing the ISCOs is much 
more favorable for the analytic solution in the case of the other 
two sequences we examined. Fig. 8 shows the errors for the evolu- 
tionary sequence that terminates at the maximum-mass nonrotating 
model. In this case, a corotating ISCO exists for some models for 
which the analytic solution is valid. The error made with the ana- 
lytic solution is 5% for the fastest rotating model, and it is con- 
sistently better than the error of the Shibata-Sasaki formula. For 
counter-rotating orbits the error for the analytic solution is some- 
what smaller than for corotating orbits (which is expected, as the 
ISCO for corotating orbits is normally at larger radii). However, 
the error is consistently larger than the error of the Shibata-Sasaki 
formula. The corresponding errors for the supramassive sequence, 
shown in Fig. 9, are very similar to the errors for the sequence in 
Fig. 8. 



5.4 Comparison to other matching conditions 

Manko et al. (2000a) also use the quadrupole moment Q in order 
to match the analytic solution to a numerical one. However, they 
redefine the parameter b as 



AR I sco/ R ,s 
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b = ±\Ja 2 + 2aMA- M 2 , 



(51) 



where A is a new parameter, with the motivation that now A mea- 
sures the departure of the analytic solution from the Kerr metric. 
We find that the above redefinition is not necessary, as it does not 
change the solution: in other words, a solution with a given b has 
a corresponding value of A. In Manko et al. (2000a) it is not men- 
tioned that Q can be set to the numerical value only for a limited 
range of the parameter b (or, equivalently, of A). Moreover, we find 
that the illustrative solutions they give do not correspond to the 
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mass corresponding to a nonrotating model of maximum mass. 



negative branch of solutions for b (that, as we have seen, are those 
relevant for rapidly rotating stars) but rather to the positive branch, 
which is closer in behavior to the Kerr metric. 

Stute & Camenzind (2002) fix the the third parameter in the 
analytic solution, b, by matching the value of the metric function 
g tt at the stellar equator. However, this is a local quantity, so that the 
analytic and numerical metrics are matched only at a single point 
in the (p,z) plane. Experimenting with this choice, we found that 
one can obtain numerical solutions that are also limited to rapidly 
rotating stars. However, since the parameter b is not fixed directly, 
but only indirectly, one has to solve a nonlinear equation in order to 
obtain b for a given value of g tt at the equator. This procedure could 
lead to multiple solutions, and one has to choose the one closest to 
a rapidly rotating neutron star by examining other properties of the 
spacetime (e.g. the higher multipole moments). As we have seen, 
fixing Q also leads to multiple solutions, however, we find that the 
procedure of fixing directly three leading multipole moments (M, 
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j and and selecting the desired solution according to the value 
of a fourth multipole moment (53) is more intuitive that fixing M, j 
and the value of a metric function a single point in the spacetime. 

In the matching procedure used by Stute & Camenzind the 
parameter b is not chosen to be real. They rather impose that b 
continuously reduces to the Schwarzschild value, b = iM, in the 
nonrotating limit (a — > 0). However, the Schwarzschild and Kerr 
solutions can only (formally) be obtained as limiting cases of the 
Manko solution by analytic continuation in the complex-fe plane. 
In a sense, the black hole solutions are "isolated points" on the 
pure-imaginary axis of the complex-fc plane, while solutions repre- 
senting neutron star exteriors lie on the pure-real axis. Therefore, 
the requirement imposed by Stute & Camenzind violates one of the 
original requirements of the analytic solution (namely, that all three 
parameters of the solution must be real). If one follows this proce- 
dure, the resulting metric components are, in general, complex. For 
complex values of b one could, in principle, use the real parts of 
some quantities in order to compute an estimate for the location of 
the ISCO, as was done by Stute & Camenzind. However, in such 
cases, even the coordinates in which the metric is expressed become 
complex numbers. Furthermore, additional multipole moments ap- 
pear that are not present in the numerical solution, rendering the 
analytic solution inappropriate for describing the physical proper- 
ties of a rotating neutron star. 

Finally, an important point in the matching procedure is to use 
the correct correspondence between the coordinates in the analytic 
exterior spacetime (10) and the numerical spacetime (1). Stute & 
Camenzind transformed the analytic metric to Boyer-Lindquist like 
coordinates, but these are not the coordinates used in (1). This can 
easily be seen when one considers that the metric (1) reduces to 
the Schwarzschild metric in isotropic coordinates (not in the usual 
Schwarzschild coordinates) in the nonrotating limit. 



6 CONCLUSIONS 

We have investigated the properties of a closed-form analytic so- 
lution for the exterior spacetime of rapidly rotating neutron stars. 
We matched it to highly-accurate numerical solutions, imposing 
that the quadrupole moment of the numerical and analytic space- 
times be the same. For the analytic solution we considered, such a 



matching condition can be satisfied only for very rapidly rotating 
stars. We found that solutions belong to two branches, only one of 
which is a good approximation to the exterior of rapidly rotating 
neutron star spacetimes. In order to evaluate the accuracy of the 
analytic solution in describing rapidly rotating neutron stars, we 
presented a comparison of the radii of ISCOs obtained with a) the 
analytic solution, b) the Kerr metric, c) an analytic series-expansion 
derived by Shibata & Sasaki and d) a highly-accurate numerical 
code. In most cases we found that the analytic solution has an ac- 
curacy consistently better than the Shibata-Sasaki expansion up to 
0(v' 4 ), for corotating orbits. Only for counterrotating orbits does 
the higher-order Shibata-Sasaki expansion perform better than the 
analytic solution. We have only shown direct comparisons for three 
constant rest-mass sequences and one representative EOS (FPS); 
however our qualitative conclusions also hold for other EOSs. 

The analytic solution we studied in this paper could become 
useful in constructing outgoing-wave boundary conditions for sim- 
ulations of pulsating relativistic stars, and for the computation of 
quasinormal modes of oscillation as an eigenvalue problem (a long- 
standing problem in relativistic astrophysics). Another potential 
application is the study of high-frequency variability in accretion 
disks around rapidly rotating relativistic stars. We emphasize, how- 
ever, that this analytic solution is only valid for rapidly rotating 
stars, contrary to previous claims in the literature. For stars of in- 
termediate rotation rates one can use the exterior analytic solution 
by Hartle & Thorne (1968), valid to second order in the rotation 
rate. This approximate solution is determined by the three multi- 
pole moments M, j and Q, but higher-order multipole moments are 
ignored. It would be interesting to determine whether the region 
in which the second-order Hartle-Thorne metric is valid to some 
accuracy, overlaps with the region in which the analytic solution 
considered here is valid. Such a study, along with a characteriza- 
tion of the spacetimes using invariant quantities (constructed in the 
Newman-Penrose formalism) will be reported elsewhere (Berti et 
al, in preparation). 
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-0.5274 


2.0800 


10.42 
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0.401 
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-14.09 


4.462 


-40.03 


-0.5443 



Table 1. Equilibrium properties for three sequences of constant rest mass 
Mb, constructed with EOS A. We show a sequence that corresponds to a 
gravitational mass of M = 1 AM e in the nonrotating limit (Mb = 1 .589M ), 
a sequence that terminates at the maximum-mass nonrotating model in the 
nonrotating limit (M B = 1.948M Q ) and a supramassive sequence (M s = 
2.O38M ). 
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Table 2. Same as Table 1, but for EOS AU. 
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-27.10 


-0.7381 


1.0578 


6.831 


1.577 


1.425 


0.0881 


13.48 




4.847 


-13.79 


2.668 


-35.06 


-0.6986 


1.0364 


7.007 


1.624 


1.428 


0.0964 


14.02 




4.699 


-15.40 


2.819 


-41.46 


-0.6629 


1.0157 


7.150 


1.672 


1.430 


0.1044 


15.05 




4.049 


-17.02 


2.961 


-48.26 


-0.6244 










M B 


2.1051/. 












3.3900 


0.000 




1.802 


0.0000 


9.276 


6.674 


6.674 


0.000 


0.000 


0.000 




2.7939 


3.204 


1.462 


1.806 


0.0073 


9.708 


4.911 


7.741 


-0.830 


1.160 


-0.556 


- 


2.5016 


4.833 


1.543 


1.813 


0.0184 


10.05 


3.904 


8.350 


-2.200 


1.847 


-2.402 


- 


2.2399 


6.184 


1.645 


1.821 


0.0337 


10.49 


2.961 


8.900 


-4.303 


2.519 


-6.549 


- 


2.0056 


7.260 


1.772 


1.831 


0.0525 


11.03 


2.087 


9.374 


-7.222 


3.185 


-14.18 


- 


1.8461 


7.887 


1.885 


1.840 


0.0685 


11.53 


1.477 


9.659 


-10.01 


3.683 


-23.05 


-0.3792 


1.6992 


8.370 


2.017 


1.849 


0.0856 


12.16 


0.881 


9.846 


-13.37 


4.181 


-35.48 


-0.5965 


1.6079 


8.613 


2.117 


1.855 


0.0976 


12.70 


0.449 


9.867 


-15.97 


4.515 


-46.22 


-0.6566 


1.5426 


8.760 


2.199 


1.860 


0.1069 


13.24 


0.035 


9.759 


-18.16 


4.771 


-55.96 


-0.6736 


1.5000 


8.846 


2.260 


1.864 


0.1135 


13.82 




9.492 


-19.80 


4.951 


-63.67 


-0.6739 










M B 


= 2.226M 












3.2103 


8.452 


1.628 


1.914 


0.0492 


9.977 


3.346 


11.04 


-5.570 


3.409 


-9.594 




3.0701 


8.350 


1.652 


1.914 


0.0496 


10.07 


3.268 


10.96 


-5.733 


3.417 


-10.03 




2.9361 


8.301 


1.679 


1.915 


0.0507 


10.18 


3.163 


10.92 


-5.991 


3.452 


-10.72 




2.8079 


8.284 


1.709 


1.916 


0.0524 


10.30 


3.035 


10.88 


-6.328 


3.506 


-11.64 




2.6854 


8.326 


1.742 


1.918 


0.0550 


10.44 


2.875 


10.88 


-6.802 


3.593 


-12.98 




2.3488 


8.607 


1.864 


1.925 


0.0665 


10.94 


2.279 


10.96 


-8.917 


3.973 


-19.46 




2.0544 


8.973 


2.024 


1.936 


0.0833 


11.62 


1.554 


11.08 


-12.24 


4.498 


-31.21 


-0.3080 


1.8375 


9.239 


2.193 


1.946 


0.1007 


12.45 


0.831 


11.07 


-16.13 


5.018 


-47.06 


-0.5536 


1.7185 


9.377 


2.318 


1.954 


0.1131 


13.27 


0.162 


10.83 


-19.18 


5.383 


-60.96 


-0.6168 


1.6800 


9.333 


2.346 


1.955 


0.1145 


13.56 




10.60 


-19.76 


5.423 


-63.67 


-0.6293 



Table 3. Same as Table 1, but for EOS FPS. 
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Q. 


/ 


M 


T/W 


R e 


h + 




M 2 


J 


S3 


b 


10 g cm 


10 s 


1 r>45 2 
10 g cm 


M Q 




km 


km 


km 


i 3 
km 


km 


1 4 
km 


km 










M B 


= 1.510M, 














0.4326 


0.000 




1.402 


0.0000 


14.83 






0.000 


0.0000 


0.000 




0.4266 


1.816 


2.162 


1.404 


0.0124 


15.14 


- 


- 


-3.671 


0.9726 


-3.554 


-0.5958 


0.4188 


2.690 


2.253 


1.407 


0.0286 


15.58 


- 


0.541 


-8.464 


1.501 


-12.65 


-0.9422 


0.4111 


3.268 


2.346 


1.410 


0.0443 


16.06 


- 


1.283 


-13.20 


1.899 


-24.96 


-0.9455 


0.4045 


3.648 


2.431 


1.412 


0.0577 


16.52 


- 


1.749 


-17.29 


2.196 


-37.84 


-0.8901 


0.3980 


3.954 


2.519 


1.415 


0.0709 


17.05 


- 


2.077 


-21.42 


2.466 


-52.65 


-0.8151 


0.3916 


4.204 


2.611 


1.417 


0.0839 


17.68 


- 


2.250 


-25.60 


2.719 


-69.47 


-0.7311 


0.3853 


4.407 


2.707 


1.420 


0.0965 


18.49 


- 


2.195 


-29.80 


2.955 


-88.00 


-0.6444 


0.3800 


4.552 


2.793 


I All 


0.1070 


19.61 


- 


1.693 


-33.46 


3.148 


-105.4 


-0.5693 


0.3790 


4.573 


2.807 


1.422 


0.1087 


19.97 


- 


1.430 


-34.06 


3.179 


-108.4 


-0.5569 










M B 


= 3.232M, 














1.4700 


0.000 




2.713 


0.0000 


13.71 


10.31 


10.31 


0.000 


0.0000 


0.000 




1.2010 


2.063 


5.009 


2.720 


0.0069 


14.20 


7.868 


12.01 


-2.656 


2.559 


-2.558 


- 


1.0639 


3.339 


5.273 


2.731 


0.0200 


14.65 


6.260 


13.22 


-8.015 


4.361 


-13.43 


— 


0.9552 


4.343 


5.591 


2.746 


0.0372 


15.19 


4.888 


14.27 


-15.78 


6.014 


-37.00 


- 


0.8692 


5.088 


5.948 


2.762 


0.0564 


15.80 


3.767 


15.14 


-25.24 


7.496 


-74.63 


- 


0.8017 


5.608 


6.320 


2.778 


0.0752 


16.45 


2.889 


15.80 


-35.49 


8.777 


-124.1 


-0.5411 


0.7495 


5.964 


6.692 


2.793 


0.0927 


17.14 


2.189 


16.28 


-45.99 


9.885 


-182.6 


-0.8287 


0.7101 


6.192 


7.037 


2.806 


0.1077 


17.82 


1.623 


16.56 


-55.82 


10.79 


-243.8 


-0.9222 


0.6729 


6.371 


7.431 


2.820 


0.1234 


18.73 


0.958 


16.65 


-67.15 


11.73 


-321.3 


-0.9373 


0.6600 


6.423 


7.586 


2.825 


0.1293 


19.19 


0.628 


16.57 


-71.63 


12.07 


-353.8 


-0.9275 










M B 


= 3.470M, 


s 












1.3847 


6.095 


5.855 


2.929 


0.0607 


14.86 


5.000 


17.96 


-24.16 


8.839 


-69.25 




1.3199 


6.030 


5.927 


2.929 


0.0610 


14.97 


4.927 


17.88 


-24.68 


8.851 


-71.55 




1.2500 


5.994 


6.021 


2.931 


0.0623 


15.12 


4.793 


17.83 


-25.72 


8.938 


-76.10 




1.1992 


5.993 


6.103 


2.932 


0.0640 


15.24 


4.657 


17.83 


-26.85 


9.058 


-81.19 




1.1160 


6.042 


6.268 


2.937 


0.0685 


15.50 


4.349 


17.89 


-29.62 


9.381 


-94.07 




0.9665 


6.282 


6.712 


2.952 


0.0835 


16.22 


3.499 


18.21 


-38.85 


10.44 


-141.1 


0.1274 


0.8781 


6.493 


7.115 


2.967 


0.0981 


16.90 


2.791 


18.51 


-48.44 


11.45 


-196.2 


-0.5157 


0.8172 


6.639 


7.482 


2.980 


0.1111 


17.56 


2.199 


18.70 


-57.65 


12.30 


-254.4 


-0.7332 


0.7604 


6.764 


7.927 


2.995 


0.1261 


18.49 


1.463 


18.75 


-69.21 


13.28 


-334.2 


-0.8412 


0.7580 


6.769 


7.948 


2.996 


0.1268 


18.54 


1.423 


18.74 


-69.78 


13.32 


-338.4 


-0.8437 



Table 4. Same as Table 1, but for EOS L. 



Approximate Matching of Analytic and Numerical Solutions for Rapidly Rotating Neutron Stars 





Q. 


/ 


M 


T/W 


R e 


h+ 




M 2 


J 


S3 


b 


10 g cm 


10 s 


1 r>45 2 
10 g cm 


Ms 




km 


km 


km 


i 3 
km 


km 


1 4 
km 


km 










M B 


= 1.551M, 














0.9950 


0.000 




1.403 


0.0000 


11.55 


0.874 


8.737 


0.000 


0.0000 


0.000 




0.9800 


2.756 


1.368 


1.407 


0.0124 


11.81 


- 


2.425 


-1.965 


0.9349 


-1.750 


- 


0.9700 


3.481 


1.396 


1.409 


0.0202 


11.99 


- 


2.856 


-3.203 


1.203 


-3.673 


-0.3370 


0.9600 


4.045 


1.425 


1.411 


0.0280 


12.18 


- 


3.195 


-4.464 


1.428 


-6.064 


-0.6044 


0.9500 


4.506 


1.455 


1.412 


0.0356 


12.37 


- 


3.470 


-5.723 


1.624 


-8.852 


-0.7155 


0.9400 


4.896 


1.486 


1.414 


0.0432 


12.58 


- 


3.699 


-7.002 


1.802 


-12.03 


-0.7689 


0.9200 


5.523 


1.551 


1.418 


0.0581 


13.05 


- 


4.036 


-9.608 


2.121 


-19.50 


-0.7935 


0.9000 


6.007 


1.621 


1.422 


0.0727 


13.62 


- 


4.209 


-12.31 


2.411 


-28.49 


-0.7644 


0.8800 


6.383 


1.696 


1.426 


0.0870 


14.41 


- 


4.129 


-15.14 


2.681 


-39.17 


-0.7089 


0.8700 


6.541 


1.737 


1.428 


0.0943 


15.04 


- 


3.855 


-16.64 


2.814 


-45.32 


-0.6744 










M B 


= 2.672M; 


3 












2.6000 


0.000 




2.205 


0.0000 


10.12 


9.404 


9.404 


0.000 


0.000 


0.000 




2.4500 


2.820 


2.298 


2.212 


0.0059 


10.28 


7.707 


10.92 


-0.905 


1.605 


-0.526 


- 


2.3000 


3.828 


2.352 


2.217 


0.0115 


10.45 


6.983 


11.42 


-1.846 


2.230 


-1.606 


— 


2.1000 


5.360 


2.451 


2.228 


0.0244 


10.77 


5.835 


12.24 


-4.059 


3.254 


-5.299 


- 


1.9000 


6.791 


2.593 


2.243 


0.0432 


11.21 


4.628 


13.06 


-7.604 


4.362 


-13.71 


- 


1.8000 


7.465 


2.687 


2.253 


0.0555 


11.50 


3.990 


13.48 


-10.06 


4.967 


-20.93 


- 


1.7000 


8.066 


2.799 


2.263 


0.0693 


11.86 


3.343 


13.86 


-13.07 


5.592 


-31.10 


- 


1.6000 


8.605 


2.939 


2.276 


0.0853 


12.33 


2.660 


14.19 


-16.82 


6.263 


-45.57 


- 


1.5000 


9.048 


3.113 


2.290 


0.1032 


12.97 


1.909 


14.42 


-21.52 


6.977 


-66.31 


-0.2276 


1.4000 


9.384 


3.337 


2.307 


0.1233 


14.17 


0.767 


14.21 


-27.54 


7.756 


-96.55 


-0.5089 










M B 


= 2.800M : 


3 












2.5000 


8.310 


2.612 


2.335 


0.0540 


10.75 


4.968 


15.03 


-9.015 


5.377 


-16.70 




2.4000 


8.310 


2.647 


2.336 


0.0556 


10.86 


4.841 


15.02 


-9.418 


5.447 


-17.93 




2.3000 


8.342 


2.686 


2.338 


0.0579 


10.99 


4.679 


15.01 


-9.976 


5.549 


-19.68 




2.2000 


8.457 


2.734 


2.342 


0.0617 


11.13 


4.450 


15.07 


-10.82 


5.726 


-22.32 




2.1000 


8.596 


2.791 


2.346 


0.0664 


11.30 


4.185 


15.15 


-11.91 


5.942 


-25.97 




2.0000 


8.799 


2.860 


2.351 


0.0728 


11.52 


3.850 


15.27 


-13.39 


6.233 


-31.09 




1.9000 


9.023 


2.943 


2.358 


0.0807 


11.78 


3.467 


15.40 


-15.26 


6.577 


-38.03 




1.8000 


9.258 


3.044 


2.367 


0.0901 


12.10 


3.030 


15.55 


-17.62 


6.979 


-47.47 




1.7000 


9.494 


3.168 


2.377 


0.1013 


12.52 


2.521 


15.68 


-20.64 


7.449 


-60.54 




1.6000 


9.705 


3.321 


2.388 


0.1143 


13.11 


1.894 


15.72 


-24.45 


7.983 


-78.60 


0.0123 



Table 5. Same as Table 1, but for EOS APRb. 
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EOS(sequence) 


Jcrit 


Q-cril /Q-Kepler 


A(1.44) 


0.39 


0.72 


A(MM) 


0.50 


0.89 


FPS(1.44) 


0.30 


0.57 


FPS(MM) 


0.50 


0.89 


L(1.44) 


0.23 


0.40 


L(MM) 


0.52 


0.87 



Table 6. Minimum (critical) rotation parameter j cr j, and corresponding ratio 
of critical angular velocity Q. crit to Keplerian angular velocity Q. Kepl „, for 
which the matching condition (36) has a real solution. 



EOS(sequence) 


AScrit 


ASKepler 


A(1.44) 


-11% 


-4% 


A(MM) 


-5% 


-2% 


A(SM) 


-2% 


-4% 


AU(1.44) 


-11% 


-5% 


AU(MM) 


-9% 


-7% 


AU(SM) 


-11% 


-10% 


FPS(1.44) 


-11% 


-4% 


FPS(MM) 


-3% 


-2% 


FPS(SM) 


-7% 


-4% 


L(1.44) 


-45% 


-4% 


L(MM) 


-3% 


-2% 


L(SM) 


-12% 


-4% 


APRb(1.4) 


-13% 


-6% 


APRb(MM) 


-8% 


-6% 


APRb(SM) 


-11% 


-11% 



Table 7. Relative difference in 53 between the negative branch of analytic 
solutions and the numerical solution, for different EOSs and evolutionary 
sequences. The difference is tabulated for the minimum (critical) rotation 
rate for which a real analytic solution exists (AS cr n) and for the model at the 
mass-shedding limit (ASxepler)- 



